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ABSTRACT 


The finite element method is applied to the 
two-dimensional, inviscid, compressible radial 
equilibriun equation for axial compressors. 
Isoparametric elements are used along with three-point 
Gaussian integration for stiffness matrix evaluation. 
The radial equilibrium equation is put into 
quasi-harmonic form for stream function formulation 
and results are presented using an isentropic flow 
assumption. Axial velocity profiles at rotor and 
stator blade edges are compared with published 
performance data of tke NASA Task-1 stage transonic 
compressor and with numerical finite element results 
of Hirsch and Warzee. 
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I. INTRODUCTION 


A. PROBLEM STATEMENT AND OBJECTIVE 


The prediction of meridional flows within turbomachines, 
be they compressors or turbines, is a difficult but 
important part of the design process. The difficulty arises 
from the presence of three- dimensional and viscous effects 
within all turbomachines and the importance arises from the 
necessity to design accurately and efficiently. 


To simplify the problem of viscous, three-dimensional 
amalysis, Wu [Re£.1] showed that this complicated flow may 
be analyzed by solving two interrelated flows: on2 being the 
blade-to-blade flow describing the flow between rotating 
blades and the other being the meridional through flow which 
describes the radial equilibrium. These flows are depictad 
in Fig 1. In addition, an inviscid and axisymmetric 
assumption is made in the through-flow thereby simplifying 
the flow to a two-dimensional, axi syametric, inviscid, and 


compressible analysis. 


Three methods may be found in current reports regarding 
the solution of the radial equilibrium equation. The first 
two are the streamline curvature method [Ref.2,3,and 4] and 
the matrix method [Ref.5 and 6] which is basically a finite 
difference technique. The third, a relatively new method, is 
the finite element method. As shown by Hirsch and Warzee 
(Ref.7]°, the solution of the radial equilibrium equation by 


the finite element method is achieved by arranging the 


equation for the stream function in quasi-harmonic fora. 


Due to the excellent results reported in R2f.7 and to 
further the research effort for finite element techniques in 
fluid flow problems, the purpose of this thesis is two fold. 
Firstly, the goal was to formulate a computer program for 
solution of the radial equilibrium equation paralleling the 
steps as presented by Hirsch and Warzee. Secondly, after 
Suitable verification of computer results with those of 
Hirsch and Warzee, the goal was to compare computer 
predicted flows with measured performance data of the Naval 


Post Graduate School's transonic compressor, 


The purpose of this paper is to present a report on the 
results obtained thus far. In Section II, the derivation of 
the radial equilibrium equation is presented followed by the 
application of the finite element method to this equation. 
Section III describes the computer program in some detail. 
Section IV contains selected test cases which wer2 used in 
program testing and checking. In Section V, conclusions are 
presented along with recommendations for further study and 
work on the project. The appendices contain the prograa 
listing along with a sample test case for reference py the 
user. In addition, a list of references is contained for 


further reading on the subject of this paper. 


Figure 1 - 


MERIDIONAL 
PLANE 


BLADE=TO- BLADE 
PLANE 


MERIDIONAL AND BLADE-TO-BLADE PLANES 


A. THE DERIVATION OF THE RADIAL EQUILIBRIUM EQUATION 


The following discussion is taken from Ref. 7 with 
slight changes in notation. The basic turbomachine geometry 
to be analyzed is depicted in Pig 2, Although the machine 
noted is one stage of a compressor, a similar analysis to 
the one that follows may be applied to other machines such. 


as axial turbines and mixed-flow machines. 


One begins with the Euler equation assuming the viscous 
forces to be negligible. 


* +(V.v)V= vPlf (II.A. 1) 


The continuity equation, assuming unsteady flow is, 


it +9 (pv) =O {IT. A. 2) 


The First Law of thermodynamics in a fluid field 
becomes, 


TY¥s = vh- vPlP (tT. A. 3 
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¢ OUTER CASING 


STATOR 


AXIS OF SYMMETRY 


Figure 2 = TURBOMACHINE GEOMETRY 


Substituting equation (II.A.3) into equation (II.A. 1) 
leads to the Crocco equation, 


~? 
st - Vs (vx) = Tys- VH (I. A. 4) 


where 4 is the total enthalpy. 


Assuming a steady and adiabatic flow, the energy 


equation becomes simply, 


(V-y)H=0O (II.A.5) 


which shows that along a streamline in a stationary systea, 
the total enthalpy is constant. 


In a relative systea, such as the case in a rotor blade 


_* 
row, the total relative velocity,W, can be expressed in the 


following form, 


—_ . . . 
where W is the constant angular velocity and J is the 


constant peripheral speed of the relative system. 


Now, the Crocco equation in a relative system becomes, 


= “a 2 2g? 
WY (ye) = Tys-Vn+ ¥ - SE) ran 


Parallel to equation (II.A.5) for the stationary system, the 
energy equation, assuming steady and adiabatic (relative) 
flow in a relative system, becomes 


‘ 
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ial 


a. (II. A. 3) 


where His the relative total enthalpy expresed as follows, 


Hp=h+ AV x uy? R2 (II.A.9) 


- 5 


From the following velocity diagraa, 


Vn > Wm 


Le 


equation (II.A.9) may be arranged as follows. 


Since, 


2 Pa om (II.A. 10) 
Wim t Wig = W* = Ve +(Vo -U4) 


then, 
ae ¥.. tp -20y, ak” (II.A.11) 

and, 
“Wie = V+ ‘Gs 2UVe (II.A.12) 
Substituting equation (II.A.12) into equation (II. A. 9) 


leads to the following relation, 


He=h ~ ¥ - Up = H- UVe (II.A. 13) 


Equation (II.A. 8) shows that He is constant along a 


streamline in a relative systenm. 


Upon integrating equation (II.A.8) between tne rotor 
inlet and outlet, the Euler equation for turbonachines is 
found, 


(IIT.A. 14) 
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It may be shown [Ref.9] that by circumferentially 


averaging equation (II.A.1), and under the axi symmetric 
flow assumption the following relation is valid, 


-Vx (vx )= T- 9S -VH +F, + Fy (II.A.15) 


where Py is the body force of the blades acting on the fluid 
and all variables are mean values along the direction of the 
circumference. Hence equation (II.A.15) is an approximation 
for axi symmetric flow. As a final note on equation 
(II.A.15), since the viscous forces were neglected in 
equation (II.A.1), there must be a force introducing the 
entropy variations along the blade. This force is 
Droportional to the pressure loss coefficient and is labeled 
Fy, the dissipative force. Pj; produces work which in turn 
produces entropy production radially along the blade. Under 
the axi symmetric assumption, entropy varies axially and 
radially only and is assumed to be proportional to the 
pressure loss coefficients [Ref. 2 and 8]. 


- Due to boundary conditions imposed on the problem and 
the axi symmetric assumption, cylindrical 
coordinates, (r,6@,z),will be used in all subsequent analysis. 
Therefore, equation (II.A.15), in cylindrical coordinates 
and with axial symmetry is as follows, 


Ve 2. i oH : 
25 (tv) ACA * +\) « eae be tte Fe ar (TE he 16} 


Ye & (ve) * ¥ > (eve) = Fe (II.A.17) 


+ > T 
Ve sae ~ 5e\4) - 2% (ey) = su. \ $$ -F, (II.A. 18) 


It is important to note here that under the axi 
Symmetric assumptions, equation (II.A.15) reduces to the 
following, 


—_ — 
V: ri? O (II.A. 19) 


Likewise in a relative system (rotor), the axi symmetric 
assumption leads to the following, 


—_ 


bed 
WwW: Pees (IT.A.20) 


Equation (II.A.16) describes the meridional through flow 
radial equilibrium equation for the finite element method. 
Since one is concerned with the meridional plane, the 


following derivative expression is taken from Fig 3. 


AX vst a. (IT.A.21) 
Von Sim = Veoe * Vo a 
Therefore equation (II.A.17) reduces to, 
. d 
RF, F Vim Sin (Rvp) (Zi oA 6 22} 


which reveals that in a duct where there are no blades and 


therefore no blade forces, angular momentum is constant 
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d along a streamline. In that case, 


As shown in Ref.9, the circumferentially averaged 
continuity equation is the following, 


Fe(4ebve) + 5; (ye b\e) 20 (IT.A.24) 


where b is the blockage factor defined by Hirsch and Warzee 
as the tangential area reduction due to the thickness of the 
blade. 


¢ 
be : ae 


(II.A.25) 


where t is blade thickness and s is blade spacing. 
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Figure 3 - MERIDIONAL PLANE 
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One further step in the formulation of the radial 
equilibrium equation for solution by the finite element 
method involves introducing the stream function. In 
cylindrical coordinates, the stream functions are defined as 


follows, 


4 3@ 
\g = hb iy (II.A.26) 
ia (i | (II.A.27) 
= yRb ot 


Substituting these expressions into equation (II.A. 16), 


=— 


the equation becomes, \ es 


Gee) + (1 st\_ 1 at +ds 
se 06 Se + 5 Gan3e |= Git -T oe (II.a.28) 


The right hand side of equation (II.A.28) is applicable 
to the absolute flows in the stator and duct regions. For 
relative flows such as those in the rotor, the right hand 
side is modified by replacing the total enthalpy,H, by the 


relative total enthalpy,H and the quantity Voe/R is 


Rr’ 
replaced by We/R. 


As a last assumption in the formulation of the governing 
relation for the meridional through flow radial equilibrium 
eguation, both the radial component of the body force, Ft 
and the radial component of the dissipative force,F,, are 
neglected. This assumption, [{Ref.1,8] does not hamper the 
accuracy of the results for conditions at design speed. 
Even though published compressor performance data used for 


19 


the test case in this thesis was obtained at 0.5 design 
speed, these force terms were also neglected in the computer 
program. As will be shown later, this assumption could 
possibly have had adverse effects on the predicted axial 
velocity profiles at the rotor hub and tip regions. 


The final representation of the meridional radial 
equilibrium equation to b2 solved by the finite element 


method is as follows, 


a (a a {2% 
(eS + (ky ~ { =O (7.4.29) 
where, 
K= ~&b (II.A.30) 


and 


erat 
YE BR Gw)] erro 


B. THE FINITE ELEMENT METHOD APPLIED TO THE RADIAL 
EQUILIBRIUM EQUATION 


In order to formulate equations (II.A.29) through 
(II.A.31) in matrix form for solution by the finite element 
method, one must apply a weight2d residual technique to the 
equations for numerical solution. The weighted residual 
method used here is the Galerkin's Method. The following 
discussion is taken from Ref. 7 with only slight changes in 


notation. 


Rewriting equation (II.A.29) and dividing through by R, 
one has, 
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hsiae (77.3. 1) 
2 (ys¥) 3 [pay = 

Kiel se)+ &KkH)-4}=0 
where this equation represents the flow in the volume,V. 


The boundary condition for this partial differential 
equation, after dividing through by &, is, 


rata # oh, (¥-¥,) =O (II.B. 2) 


where this equatioa solves the flow on the closed boundary 


of the volume,or,sS. 


By applying the weighted residual process to equations 
(II.B.1) and (II.B.2) and using an arbitrary weighting 


function, W(r,Z), one has 
\w (G2) tone + [wit2) tue dS =O (TT. Bs 3} 
V s 


where See and D coy are the volume and surface residuals 


respectively, or, 


i | 
, EF st) + (x) . f 3° a atari dg 


Y 
Ys * rates + ff, (Y- 4.) (II.B.5) 


If the solution to equation (II.B.1) was exact, both Fs 


and cr would be equal to zero. 


suc 


In order to clarif the boundary condition, equation 
y 


(II.5.2), one may analyze the equation as follows. 
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On the surface,S,where is specified, 


Y- y (II.B.6) 


and, 


4, -* co (II. B. 7) 


oY 
Similarly, on the surface, where yi ?,s where 


a 


(II.8.3) 


&. 
l 
1e' 


mt Se =O, S, US, “.% (II. B. 9) 


Due to the axi symmetric assumption, the final equation 
will not involve dv and dS but the intersection of d¥ and 4s 
with the meridional plane. Therefore, one must transform the 
volume integral,dV, to a surface integral and the surface 
integral,dS, to a line integral. 


Hence, let, 
dr 
dC 


intersection of dV and meridional plane 


intersection of dS and meridional plane 


mae dN= 2ur Rd 
dS= Iredc 


With this transformation, one may rewrit2 equation 
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(II.B.3) as follows, 
(-woea) |e) « 3k a \s {|2n dQ + 
n weak Inde =O 


where on the contour, og Y=. 


(II.B.10) 


One must now integrate the first term in equation 
(II.B.10) by parts to obtain the following, 


(wie ke) ky Venda. res 


av y 
+ (k oie Be 3 J + Jed ‘t*MdC=O 


Inspecting the first term in eguation (II.B.11), one may 


(II .B.11) 


use the following integral theorem to simplify further 


[Lagan = (onde (II.B. 12) 
Cc 


we 


Rewriting note (ZE.B.11) ,gives;, 


Vy, av o¥ IW 
- hu e" 3m oC- fd me? ot 3) ae 
foe andl = O 
Pinally, since 
ca = $F ite + 353" (II.B. 14) 


equation (II.B.13) reduces to the following, 
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ase BM )-{W]da=o (II.B. 15) 


One now has the final equation in the form for use by 
the weighted residual method using any arbitrary weighting 
function, WW (r,2). As noted previously, the Galerkin's 
Method will be used here which implies that the weighting 
functions are the same functions used in approximating the 


stream function, ¥ ; 


Before applying the finite element method, one nust 
discretize the continuum and then approximate the unknown 
function ¥ by a set of polynomials. For this particular 
problem, eight-noded iso-parametric elements were chosen for 
discretization, see Fig 4, and the following approximatina 


functions were used. 


g 
Y= 2 Ni (Gn)¥ (IT.3. 16) 


where, 
Ni (&) = shape functions 
Y, 


Y= value of YW at any arbitrary 


value of Y at the node 


location within the element. The shape functions,N:, used 
here are defined by the following relations as shown in 
Ref.10, 
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Ni (Em) = E(l+ GEC yelE Gs My ~ 1) 
Ni(4,y)= 2 (\- @*)( le) (II.3.17) 


Ne Gm) = E (t+ FRAC) 


where the following coordinate transformations are used, 


8 
Cs 2 NG 4 (II .B.18) 


: 
z= 2 Ni (4)¥)8 


3) 2 i! 


1 
rash 
> 
| LS 
—e 
co 
ee 
i 
+ 
fea 


uw 
fon) 
~N 


Figure 4 - ISOPARAMETRIC QUADRILATERAL ELEMENT 


At this point, one is ready to apply the Galerkin Method 
to equation (II.B.15) by substituting equation (II.3.16) for 
the unknown function Y , and N; for the weight function, W, 


which yields, 
(x fe Syl), > 4. OM) dQ - rr dn = Ot.8.19) 
wD url L3 = 


This integration yields the following system of 


equations which is solved for the unknown nodal, 


(II.B.20) 
Kwi Kw Yn Ee 
where, 
ONC ON ONS UN: 
Ki = k se +> a dz (II.B.21) 
3 
and, 


f= [EN de (II.B.22) 


In addition since both the "stiffness matrix',K, and the 
right hand side vector, [FJ], are functions of Y , the system 
as defined by equations (II.B.20) through (II.8.21) must be 
solved iteratively. 


At this point one has the total finite element 
formulation of the radial equilibrium equation as defined by 
equations (II.B.19) and (II.B.20). The problems which 
remain to be clarified are basically two fold. Firstly one 


za? 


must evaluate the integrals in equations (II.5.19) and 
(II.B.2) by numerical methods, and secondly, th2 solution 
procedure for the non-linearity must be formulated. In Part 


C, both of these final steps are presented. 


C. NUMERICAL INTEGRATION OF STIPPNESS MATRIX AND SOLUTION 
PROCEDURE 


As noted in Section II.B, e?2valuation of equation 
(II.B.21) must be perfocmed numerically. In addition, one 
realizes that the derivative expressions enclosed within the 
interval must be evaluated by a coordinate transformation. 


This is done in the following way, 


Since, 
g 
2 Nil) 
(i7.C. 1} 
$ 
t= a Ni l§.y)% 
then, 
dai | Nee, ON de 
fe r 
i t a {TiC 2 
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and in matrix fora, 


di dt dr |{ SNe 

> > 
+4 rss { 4 . (EE .CeS) 
2Ni 4 dc }] aN 


Furthermore, defining the Jacobian matrix as, 


ear 
J = “4 *4 (t1.c. 8} 
2 oar 


ee 


then by dividing both sides of equation (II.C.3) by J, one 


has the following transfocmation, 


me 3 
si S (x]" "7 (Iz.C.3) 
oN oN. 


In addition, it has been shown [{Ref.9] that 


atdr |s| d¢dy (II.C.6) 


Now, with equations (II.c.5) and (TE.C.5)» equation 


(II.B.21) becomes the following, 
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c= S(kl3 ae sj hes) by iat" | ately (II.C.7) 
Equation (IZ.c.7) is best eee using the 


Gauss-Legendre integration method since it is of the 


following form, 


Ki [Jotmegey (II.C.8) 


or finally, £ Ref. 10}, 


; 4 14.8; 5 (8, r) eye 


where A and B.are coefficients (Fig 5) for both two and 
. 


three point Gaussian Quadrature. 


At this point, one has the tools to calculate all 
the elements of the stiffness matrix. In like manner, the 
right’ hand side vector,f, is calcualted by numerical 
integration. 
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NUMBER OF 
GAUSS IAN 
POINTS 


0.57735 02691 1.00000 00000 


0.77459 66692 0.55555 55555 


0.00000 00000 0.88888 88888 


Figure 5 - GAUSSIAN INTEGRATION POINTS 
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2. Solution procedure 


The following is a synopsis of tthe basic solution 
process. Specific details concerning equations and methods 
of computer coding are covered in the proceeding section. 
The proceeding is meant to give the reader a preview of the 


solution process. 


a. Discretization 


Initially the machine under analysis is 
discretized into eight-node iso-parametric elements. The 
axial calculation stations are placed arbitrarily in the 
duct regions and along blade edges and centers for the rotor 
and stator as shown in Fig 6. At this point the  syst2n 
topology and nodal coordinates are specified. 


b. Initialization 


To begin the iteration process, one must assume 
an initial internal strean function, velocity, and density 
distribution. In the program, the initial internal stream 
function was assumed to be that of the outer boundary 
throughout while the velocity and density distribution was 


assumed to be that of the inlet. 


US 


| 
a 


hy 
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COMPRESSOR DISCRETIZATION 


Figure 6 - 
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c. Calculation of thermodynamic variables 


Before calculating the right-hand side vector,f 
, one must obtain distributions of angular momentua, 
enthalpy, and entropy. This is done by first calculating 
the thermodynamic variables at the inlet axial station from 
the given inlet conditions. In order to proceed axially 
through the machine to calculate the nodal angular momentun, 
enthalpy, and entropy, tne following three equations derived 


in Section III are used. 


H = CeT = constant along a stator streamline 
2 
(wr) 
He= Cle “5 = constant along a rotor streamline 


v Ve = constant along a duct streamline 


An example of this calculation procedure for the jiuct region 


is shown graphically in Pig 7. 
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Figure 7 - 


DUCT ELEMENT 
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In this figure, the angular momentum at point X 


is equal to the angular momentum at point Y. More formally, 


§ 
(rVg), = = NG m) {rVe). 2 (ve) y (12.0.2. 


Since the previous axial station's thermodynamic 
“variables are known, one must now find the values of § and 
n at point Y. This is done iteratively in th2 following 


way. Since 
S 
se Y| = 2 Ni (4,4) ¥: (F8.6.2.0 
u=3 
and along the left side of the element, 
&= -{ it.c2H 
{ 


then equation (II.C.2.2) may be solved for i by a suitable 
iteration method. As will be shown in the next section, a 
half-interval method was used to obtain the unknown YW. 
Once al is known, then equation (II.C.2.1) is solved for 
the angular momentum at point X. The rotor and stator are 
handled in a similar fashion. In addition, the rotor and 
stator deviate the flow creating a three-dimensional flow 
field between the blades in the respective blade row. Low 
speed cascade correlation data [{Ref.13] was used to 
calculate the effective turning angles in the rotor and 
stator. These effects are calculated beforehand with known 
mass flow rate and unifora axial velocity assumptions at the 
rotor inlet. The results of these calculations are part of 


the input data routine in the form of relative and absolute 


flow angles at the rotor nodes and absolute flow angles at 


the stator nodes. This will be shown more exactly in the 


next section. 
4d. Calculate matrices 


At this point, the right hand side vector, f, and 


the stiffness matrix ,K, are calculated. 
@. Solve system of equations 


The system of eguations as shown in equation 
g 3 


(II.B.20) is solved for the nodal stream FUnCELON. 
£. perform relaxation iteration 


Due to the strong non-linear properties of tae 
system of of equations, the following iterative scheme is 


necessary. 


gga [Fyn] arene 
& & 


t c 


where d is the under relaxation factor. As will be shown in 
Section III, tnis scheme is performed only in certain 
regions of the machine and in addition after a specified 


number of iterations. 
yj. Update velocity and density profiles 


Using the current nodal distribution of the 
stream function, axial and radial nodal velocity components 


are calculated along with a new nodal density distribution. 


v4 


At RA ersnssesstennvesitentette 


Again, this calculation procedure will be shown in the next 


section. 
h. Test for convergence of y 


Stream function convergence criteria is now 
tested and will determine if further iterations are 
necessary. The solution is said to converge if tha following 


equation holds for all nodes. 
ik ete cee Zé (II.C. 2.5) 
where € is a designated requirement for convergence. 
i. Summary 


In summary, tne eight steps involved in the 


solution are noted below; 
(1) Discretize the continuun. 


€2) Assume an initial stream function, velocity, 


and density solution. 


(3) Calculate the nodal thermodynamic variables 


from the given inlet conditions. 


(4) Form the right hand side vector, f(r,z), and 
the stiffness matrix, K. 


(5) Solve the system of equations, given by, [kK] 


= {F] for a new stream function distribution. 
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(6) Perform relaxation iteration if raqguired. 


(7) Calculate new nodal velocity and density 


distributions from the current stream function sclution. 


(8) Test the solution for convergenc2, and if 
required, repeat steps (3) through (8) using the current 


nodal stream function values. 
This concludes the solution description and now 


one is ready to more completely understand the computer 


program which assembles the preceding eight steps. 
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IIIT. THE PROGRAM 


A. OVERALL FLOWCHART AND DESCRIPTION 


The overall flowchart of the program is depicted in Fig 
oe Those blocks denoted by the letter 'S'" are subroutines, 
while the remaining calculations are an integral part of the 


main program. 


After proper dimensioning of all arrays and subsequent 
initialization, the input data are read and then printed. 
This not only presents a physical picture of the problem but 
also serves as a cross check to the user for correct data 
insertion. In addition, a subroutine is available to obtain 
a computer drawn plot of the mesh (Fig 6) and is a further 


check on proper data input. 


At this potnt ail the necessary variables have been 
stored and the iteration counter for strean function 
convergence is set. With the current nodal values of p and 
the given inlet thermodynamic conditions, the thermodynamic 
variables throughout the machine are calculated. From the 
calculated values of enthalpy and angular momentum, 
(isentropic flow is assuned), the right-hand side vector is 
calculated followed by the stiffness matrix calculation 
(equation II.B.21). 


The system 9f equations (equation (I1.8.20)) is now 


solved for the new nodal stream function distribution. [t 


is here where for all iterations but the first that a 


relaxation factor is applied as noted previously in equation 
CEE SO. 204). The reasoning behind not applying the 
relaxation scheme to the value of nodal Y after the first 
iteration is the fact that the first iteration produced a 
close approximation to the correct strean function 
distribution. With this close approximation to the stream 
function came a velocity and density distribution which ina 
turn was near the correct solution. It was found that if 
the first iteration was relaxed, the second iteration became 
unstable since in fact the velocities and densities were 
themselves farther from the true values than were assuaed 


initially. 


After testing the nodal stream function for convergence 


by use of equation (II.C.2.5), the calculation proces is 


un 


either repeated or ceased by virtue of convergence or 


limiting the number of iterations. 


As stated previously, low speed cascade correlation data 
(Ref.13] were used to calculate turning angles in the blade 
regions. These angles were assumed constant throughout the 
solution and not refined after subsequent iterations. 
Further work on the sSomputer program could entail an 
additional computational routine which would calculate the 
new turning angles after each iteration. A sample 


calculation of rotor turning angles is shown in Appendix D. 


In the following sections the program structure is 


examined in more detail. 
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CALCULATE THERMODYNAMIC 
VARIABLES FROM GIVEN 


CALCULATE RIGHT HAND 
SIDE VECTOR FROM GIVEN 
THERMODYNAMIC VARIABLES 


CALCULATE STIFFNESS MATRIX 


SOLVE SYSTEM OF EQUATIONS FOR y 


PERFORM RELAXATION 
ITERATION IF NOT lst ITERATION 


CALCULATE NEW VELOCITY AND DENSITY 
DISTRIBUTION FROM NEW STREAM FUNCTION 


STREAM FUNCTION 
CONVERGENCE SATISFIED OR 
ITERATION LIMIT 


NO REACHED YES 
INCREMENT ITERATION COUNTER PRINT FINITE 
ELEMENT RESULTS 


Figure 8 - PROGRAM FLOWCHART G roi 
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B. THE MAIN PROGRAM 


1. The input routine 


The following is a description of the input data 


required by the program. The data are arranged into 


categories described in the following manner. 


a. category 1 


Problem identification. 


' b. category 2 


Number of nodes and number of elements. 


c. category 3 


Node numbers, nodal coordinates and 


blockage factor. 


d. category 4 


System topology. 


@. category 5 
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twelve 


nodal 


function. 


Element type; duct, rotor, or stator. 


category 6 


Absolute flow angles for rotor and stator nodes. 


category 7 


Relative flow angles for rotor nodes. 


category 8 


Inlet thermodynamic quantities. 


category 9 


Physical constants for fluid under observation. 


category 10 


First estimate of internal stream function. 


category 11 


7) 
+ 
"mt 
w 
w 
=] 


Node numbers and specified nodal 


category 12 
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Node numbers where the right hand side, f(r,Z), 
is to be calculated. 


Before describing in detail the format to be 
followed for data insertion, it is important to note the 


following assumptions. 


(1) Uniform flow conditions at inlet and 
outlet. 


(2) Uniform flow conditions at rotor inlet 
for calculation of appropriate turning angles. This 
assumption is necessary to calculate the values of rotor and 


stator flow angles. 


With this in mind , the discussion will 
continue. 


The following describes each category in more 
detail. 


Category 7%: 


Pormat: (20A4) 


Number of cards: 1 | 


Procedure: Enter the title of the problem in 
columns 1-20. 


Category 2: 


Format: (2510) 


Number of cards: Egual to the number of 
nodes in the system. 
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Procedure: Enter the ‘number of nodes in 
columns 1-10, and the number of elements in 11-20. Both 


integers must be right justified. 
Category 3: 
Format: (I10,3F10.0) 


Number of cards: Equal to the number of 


nodes in the systen. 
Procedure: Each card contains the node 
Number followed by the Z coordinate, R coordinate, and nodal 


blockage factor. The coordinates are in dimensions of 


inches. 
Category 4: 
Format: (915) 


Number of cards: Equal to the number of 


elements in the system. 


n 
an 


Procedure: Zach card contains nine intege 
Fight justified in columas 5, 10, 15, etc., through 45. nh 
first integer is the elenent number followed by the eight 
nodes associated with that element. It is important to note 
that the nodes are read in starting with the upper right 
hand node and proceeding in a counterclockwise fashion 


around the element. 
Category 5: 


Format: (2110) 


Number of cards: Equal to the number of 


elements. 
Procedure: Enter the element number in 
columns 1-10, followed by the integer ‘1! (duct), *2% 


(rotor), or '3' (stator) describing the element as either in 


a duct, rotor, or stator region. 


Category 6: 


Format: (5X,A4%,110,F10.0) 


Number of cards: Equal to the number of 


rotor and stator nodes plus one 'STOP' card. 

Procedure: Enter the node number (right 
justified) in columns 11-20 followed by the value of the 
associated absolute flow angle in radians in columns 21-39. 
The last card in this category is a 'STOP' card entered in 
columns 7-10. 

Category 7: 


Pormat: (6X,A4,110,F10.0) 


Number of cards: Equal to the number of 


rotor nodes plus one ‘'ST9P' card. 

Procedure: Enter the node number (rigat 
justified) in columns 11-20 followed by the value of the 
associated relative flow angle in radians in columns) 21-390. 
The last card in this category is a 'STOP' card. 


Category 8: 


Pormats: (7710.0), (F 10.0) 
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Number of cards: 2 


Procedure: Enter the following quantities in 


the prescribed order and with the noted dimensions. 


First card 

Mass flow rate: (lbm/sec) 

Inlet axial velocity: (ft/sec) 
Outlet axial velocity: (ft/sec) 
Inlet total density: (Lbmsft® ) 

Inlet static density: (lba/tt? ) 
Inlet total pressure: (1b£/in ) 
Inlet total temperature: (°R) 


Second card 


Category 9: 


Pormat: (3F10.0) 


Number of cards: 1 


Procedure: Enter the following quantities in 


the prescribed order. 


Gas constant: (ft-lbf/lbm-°R) 
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Ratio of specific heats 


Constant pressure specific heat: 
(BTU/1Lbm-°R) 


Category 10: 
Pormat: (F10.0) 
Number of cards: 1 


Procedure: Enter the first estimate of the 


internal stream functiion to be used in the first iteration. 
Category 11: 
Format: (5X,A4,110,F10.0) 


Number of cards: Equal to the number of 
nodes having a specified value of the stream function plus a 


'STOP' card. 


Procedure: This set of cards allows the 
stream function boundary conditions to be read in. A 
typical card contains an integer, right justified in columns 
11-20 , which is the node number, followed by the value of 
the specified stream function in columns 21-30. The last 


card is a "STOP card. 
Category 12: 
Pormat: (5X,A4,110) 


Number of cards: Egual to the number of 


nodes where the right hand side is to he specified. 
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Procedure: Enter the node number, right 
justified in coluagns 11-20, where the right hand side is to 
be calculated. Again, the last card in this category is a 
‘STOP! card. 


After all the data has been read by the program, 
the input data is printed and the mesh is plotted for 
verification by the aser. The samapl2 format is shown in 


Appendix C. 


This concludes the input Co 
section describes the calculation of the 5s 
K. 


De Stiffness matrix avaluation 


DtszsesP Err SSe5=22 2-255 =- === 


As shown previously in Section I¥.€.%, the Loliowzng 
equation describes each term in the @ight by eight elemental 


matrix. 


Ki noe oN de Coy det (Tag dy (TIT. B..2. 1) 


o{ ot 


In addition, 'k' is defined in the following way in order to 


numerically integrate th2 equation. 


1 
E geNilGin): Ze NlGy) b 


where b is defined as the elemental biockage factor taken as 


K = 


(LIT Be eee 


an average over the eight nodes of the particular element 
and (qm) are the defined Gauss-Quadrature integration 


points. 
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FIND AVERAGE BLOCKAGE FACTOR 
OF THE ELEMENT, BBAR 


FIND SHAPE FUNCTIONS AT 
GAUSS-QUADRATURE POINT "I" 


FIND JACOBIAN AT "I", (JJ 


FIND DENSITY AND "R" COORDINATE AT "T' 


NEXT INTEGRATION 
POIN 


FIND DETERMINANT OF JACOBIAN AT "I'' 


FIND 3_N 


3 tc oe AT "I" 


FIND 4, arr 
WN. 
me 


FIND [y)7) ar "3" 


FIND {ryt T arovr 


PERFORM INTEGRATION AT "'L" 


II G.T. NUMBER 
QF ELEMENTS 


ASSEMBLE STIFFNESS MATRIX 
W/O REGARD FOR BOUNDARY 
CONDITIONS 


STIFFNESS 
Le = LE 

MATRIX COMPLETE 
Figure 9 - STIFFNESS MATRIX EVALUATION 
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Fig 9 depicts the flowchart for both elemental 
stiffness matrix evaluation and the assemblage int 

system stiffness matrix. More specifically, the figure 
showS a three-point Gaussian Quadrature scheme but can be 
Changed to a two-point scheme by simply integrating four 


times instead of nine as shown. 


The actual coding of the stiffness matrix evaluation 
and assemblage may be found in lines STR03510 tarough 


STROY770 in the computer progran. 


At this point, th2 system of equations ar2 modified 
for the boundary conditions and solved for the nodal strean 
function values. An equation solving routine, DSINOQ, 
available in the system library was used for this purpose. 
It was found that no comparable savings was realised by 


using a banded equation solver. 


As noted previously in Section II.C, a relaxation 


scheme is necessary for convergence to a solution. 


Two distinct dirtferences with regard to the 
iteration method were noted from that of Ref.7. Firstly, it 
was found that relaxation was necessary only in the rotor 
and stator elements and also in the duct region between the 
rotor outlet and stator inlet. Secondly, due to the extreme 
non linearity in the rotor-stator areas, a switch was 
required which changed the sign of d in equation (II.C.2.4) 


as required for stability of convergence. Clarification of 
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this change follows: [It was found that during the initial 
three or four iterations, the stream function values of the 
rotor-stator nodes sometimes exceeded the value of the upper 
boundary. Due to an absc2nce of sources within the domain 
of solution, this occurrence was incompatible with the 
boundary conditions. At this point, it was necessary to 
make & negative in equation (II.C.2.4). During subsequent 
iterations, as the solution converged, the rotor-stator 


regions became stable and the sign of & was returned to its 


a 
oO 


positive value. This iteration proved to stabilize t 
solution with respect to stream function values and 


velocities. 


The iteration procedure is coded in tae comouter 


program from lines STRO5070 through STRO5200. 


Once convergence is obtained or the number o9£f 
iterations have reached the limit imposed by the user, the 
results are displayed. A sample output is shown in the 
Appendix. In addition, the units of all dependent variables 


are the same as those noted in the input routine. 


C. THE SUBROUTINES 


The following describes each of the six subroutines in 
the computer program. Each subsection contains 2 list of 
calling arguements and for subroutines FCAL, SLINE, and VEL, 
a basic flowchart. In addition, for those subroutines whose 
mathematical theory was not presented in Section III, a 


brief treatment is also given. 


a3 


This subroutine calculates the shape functions 
(equation (II.B.17)) at the values of 4 and a as requested 


in the arguement list below. 


SUBROUTINE SHAPE (£,Z,SF) 


ie] 
" 


value of 5 (input) 
Z = value of W (input) 


SP = eight by one vector of the 2ight shape functions. 


fe") 


JACOB calculates the Jacobian matrix as define 
equation(II.c.1.4) for the value of % a | denoted int 
arguement list. 

SUBROUTINE JACOB (51,21,D,2,8CS$,Z2C$,RJAC) 

E1 = value of y (input) 

Z1 = value of § (input) 
We 

D = eight by one vector of 34 (calculated) | 
: oN 

E = eight by one vector of Aiea la man 


RC$ = eight by one vector of the 'r’ coordinates of the 


nodes associated with the element (input) 


ZC$ = eight by one vector of the 'z' coordinates of the 


nodes associated with the element (input) 


RJAC = two by two Jacobian matrix (output) 


ff 
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In addition, the subroutine assumes that the vectors 
RC$ and ZC$ contain element coordinates arranged in a 
counter clockwise fashion beginning with the upper right 


corner node. 


This subroutine calculates the thermodynamic 
variables throughout the machine given the inlet conditions 
as described in Section I[I.C.2. The calling arguements are 


defined below. 
SUBROUTTNE SLINE(UINLET,RC,PSI,WRL,H,UVEL,VVEL,TVEL, 
NODE,NNODEI,CP,(cT,KK,ALP,WG,TWEL,BE,HS) 
UINLET = Inlet axial velocity 
RC = Nodal 'r' coordinates vector 


Nodal stream function vector 


mu 

n 

Lan) 
u 


= 

ee) 

~ 
i 


Nodal angular nomentum vector 


H = Nodal total enthalpy vector 


UVEL = Nodal axial velocity vector 

VVEL = Nodal radial velocity vector 

TVEL = Nodal absolute tangential velocity vector 

NODE = Matrix containing nodes associated with the 
element 

INLET = Vector containing node numbers at inlet station 

NNODEI = Number of nodes at inlet station 

CP = Specific heat 


tT 


Total temperature at inlet 


KK = Iteration counter 

NTE = Element type vector 

ALP = Nodal absolute flow angle vector 

TWEL = ‘Seda Extaiias tangential velocity vector 
BE = Nodal relative flow angle vector 

HS = Nodal static enthalpy vector 


As shown in Fig 10, the basic calculation procedure 
begins with calculating the required energy and momentun 
values at the inlet station. At this point, beginning with 
element one, the element type is interrogated to jiistinguish 
between duct, rotor, and stator elements. If the element is 
in a duct region, then the streamline intersections for 
local nodes 2,6,7,8 and 1 (Fig 7) are determined along with 
the associated values of energy and angular momentum. For 
the rotor and stator elements, one must initially find the 
energy and momentum values at local nodes 3,4,5 (Fig 7) due 
to the discontinuities imposed by the blade adges. Once 
these calculations are performed, then the process for the 
remaining nodes in the element proceeds in a similar fashion 


to the duct elements. 


After all the elements have been cycled through, the 
new distributions of nodal angular momentum and energy are 
returned to the main program for further computations. 
Specifically, these values will be used by the next 
subroutine , FCAL, for calculation of the right hand side 


vector. 


FIND ANGULAR MOMENTUM, 
TOT ENTHALPY,STATIC ENTHALPY, 
AT_ INLET STATION 


BEGIN WITH Ist ELEMENT 


DUCT <Fire net STATOR 
ROTOR 
FIND H ,rW ,W 
e ° o 
AT LOCAL NODES 3,4,5 


FIND STREAMLINE 
INTERSECTIONS OF 


FIND H,rV_,V 
e 6 
AT LOCAL NODES 3 


FIND STREAMLINE 
INTERSECTIONS OF 


FIND STREAMLINE 
INTERSECTIONS OF 


LOCAL NODES REMAINING NODES REMAINING NODES 
2,6,7,8,1 AND AND CALCULATE AND CALCULATE 
CALCULATE H, H»tW ,W 5h H,rv5,V 59h 


rv_,AND h. 
e 


NO 


ALL ELEMENTS 
HAVE BEEN 
CYCLED THROUG 


NEXT ELEMENT 


oe! 


Figure 10 - SUBROUTINE SLINE 


57 


FCAL calculates the right hand side vector as 
defined by equations (II.A.31) and (II.B.21). Using the 
identical coordinate transformations for numerical 
integration as described in Section II.C, the final equation 
to be coded is the following, 

£ | 
t + Ne ON Ne: (Lr, 3Ne 
e ZN Vj. oN ; ( ! 25 
~1-! {TI7T.c.4. 1) 


tT (2,2): ode -HiN dT dG oy 


where isentropic rlow is assumed, and, 
We = anguiar wonentun 


(IITT.C.4..2) 


total enthalpy 


a 
u 


The arguement list is defined below. In addition, 
only those variables in the list which have not been defined 


previously are described. 


SUBROUTINE PCAL(F,W,H,ZA,EA,UVEL,RC,ZC,a#&L,TVEL,NPS 
,NODE,NN,NE,NNFSP, TWEL, NTE) 


F = Right hand side vector, f(r,z) 

W = vector of gaussian quadrature coefficients : 

ZA = Vector of 8. gaussian quadrature points 

EA = Vector of gaussian quadrature voints 

NFS = Vector containing nodes where the right hand sile 
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is to be specified 


NN = number of nodes 

NE = number of elements 

NNFSP = Number of nodes where the right hand side is 
specified 


Fig 41 depicts the basic flowchart for the 
Subroutine. To initialize the procedure, one begins with 
the first node (upper right hand corner) of the first 
element. A switch is then applied which determines if the 
right hand side is to be calculated at the node or if a 
stream function value has been specified. This information 
is transferred from the main program through th2 arguement 
list. Once the node is allowed through the switch, then the 
integration process is started at the first integration 
point. As in Section III.A.2, the flowchart depicts a 
three-point Gauss Quadrature scheme. After the integration 
has been completed, a switch determines if all the local 
nodes in the element have been cycled through and if so, 
then the assembly of the slemeatal vector, FS, is performed 
to build the system right hand side vector , F. Finally, 
the subroutine determines if all the elements have been 
examined in order to signal completion of the right hand 
Side vector. At this point, the vector, F(r,z), is returned 


to the main frogram for problem solution. 
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NEXT NODE 


ASSEMBLE 


RHS VECTOR 
SYSTEM, F 


FOR 


BEGIN WITH lst 
ELEMENT, IT = 1 


Ti = [J+t 
BEGIN WITH Ist 
NODE OF ELEMENT, I = 1 


S F(R,2) TO B ——e T= I¢l 


SPECIFIED AT NEXT NODE 
THIS NODE ? NO 


YES 


BEGIN WITH FIRST 
GAUSSIAN INTEGRATION 
POINT, 4 ,>N, J = 1 


FIND SHAPE FUNCTIONS, 
JACOBIAN, INVERSE OF 


JACOBIAN DET J 
} ee ee er 
FIDE, E Zi 2 i el 
ZNR, 
3 NEXT GAUSS 
FIND = (ue 4)W, POINT 
i=l of 
a 
FIND = (L 4)H, NO 
i=l il 
FIND ELEMENTAL 
RHS VECTOR FS(L) 


NO YES 


YES a 
f+ NEXT ELEMENT 


RETURN 


Figure 11 = SUBROUTINE FCAL 
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5. Subroutine vel 


This subroutine calculates axial and radial 
velocities and also densities at each of the nodes from a 
known stream function distribution. As noted previously in 
Section I1I.C.2, both velocity and density profiles are 
updated after obtaining the latest value of nodal stream 
function. 


The velocity calculation proceeds from the strean 
function equations, 


: 
| 


Va 


TTY .0.5. % 


Vue prb a CELT. Ce Fe 2) 


where 'b' is the tangential blockage factor. Since fr, ,and 


are of the following forn, 


+“ } fo Ni (727.0. 5.5) 


tt 
§ 
He oN 
t 
then the equation for the axial velocity, Yee becomes, 
1 3 
Ve= WW 7 \2 dN yy (III.C.5.4) 
Ba PN 2 Nev, Lies OF 
7 7 


Again, since the shape function, Be AGAR is is not an 
\ 
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ERS ea eee ea ee 


implicit function of 'r' and 'z*", one must use equation 
(II.C.1.5) to obtain the proper derivatives for computation 
of equation (III.c.5.4). Por example, from equation 
HIT CLS}, 


GNA 


3 
ae (T"'(,y 1) Bo +3702) 4 (III.c.5.5) 
eFl 


_ | 


At this point, with equation (EITC. 5.9) 
substituted into equation (III.C.5.4), one has the complete 
expression for the axial velocity as functions of 8 ae One 
proceeds similarly for expressing the radial velocity, 7 , 
in terms of § andy - 


In order to calculate the nodal density, one uses 


the following density relation for flows in the stator and 


Guct regions. 
(TIT. C256) 
where Py is the stagnation density. 


Since, 


\': (vy + Vg) (I+ tin’) (III.c. 5.7) 


then, 


_—_— 


4-(1- 7 (er) (ys 4.) (\+tu2a) i lal 
ry 
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Since the density appears on both sides of the 
equation, the new nodal density is obtained iteratively at 


the node. 


For the relative flows in the rotor, the following 
relation for static density is used [{Ref.14)}. 


= 
7 2 (en — (4 w wet | ot (IET.c.5.9 


Again, the solution of the nodal density is obtained 


in an iterative fashion. 


In the following arguement Lest, only those 
variables not defined in the previous subroutine 


descriptions are noted. 
SUBRCUTINE VEL(NE,NN,2C,NODE,G,RG,TT,RHOT,RHON,ZC, 
PSI,RHO,8,UINLET,UVEL, VVEL,RdAOSTA, NTE, ALP) 
G = Ratio of specific heats 


RG = Gas constant 


RHOT = Total density at the inlet 
RHON = Work vector which contains the new nodal density 
distribution 


RHO = Nodal static density vector 
B = nodal blockage factor vector 
RHOSTA = Static density at the inlet station 
The basic flowchart for SUBROUTINE VEL is shown in 


Pig 12. Beginning with the first node of the first element, 


the Jacobian matrix (2quation (II.C.1.4)) and its inverse 
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are found. At this point the partial derivatives with 
respect to ‘r' and 'z' of the shape functions are found as 
noted in equation (III.c.5.5). A switch then allows those 
nodes not at the inlet station to pass and calculates the 
new density and velocities at the nodes. For those nodes at 
the inlet, the velocities and static densities are retained 
at the given inlet conditions. This is done to maintain 
boundary condition integrity for the solution. After 
cycling through all elements, the subroutine returns the new 
nodal velocity and density distributions to the aa 


program. 
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BEGIN WITH lst ELEMENT, [I=1 


BEGIN WITH lst NODE OF ELEMENT 
[=1 


FIND JACOBIAN 


FIND INVERSE OF JACOBLAN 


=UINLET 
v_=0 


= one 


<> NO NEXT NODE 


YES 
<test> = NEXT ELEMENT 
YES 


Figure 12 = SUBROUTINE VEL 


6. Subroutine mplot 


This subroutine utilizes the Calcomp plotter to 


depict the mesh topology of the machine under observation. 
SUBROUTINE MPLOT(RC,2Z2C,NODE,NN,NE) 
This completes the description of the main oprogran 


and associated subroutines. In the next section, a test 


case is carried through from data input to final results. 
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The progran was tested by using published performance 
data [Ref.12] of the NASA Task-1 stage transonic compressor. 
The compressor was discretized into twenty-eight elements 
and 107 nodes with 15 axial calculation stations (Fig 6). 
The speed was 0.5 design speed with a mass flow of 107.6 
lbm/sec. In addition, uniforn flow was assumed both at the 
inlet and outlet stations. Turning angles for the rotor and 
stator were pre calculated assuming uniform conditions at 
the rotor inlet and using NASA SP-36 blade correlation data 
[Ref.13]. These absolute and relative flow angles were 
assumed constant throughout the iterative procedure a hey 
were an integral part of the input data. Th2 Appendix 
contains a listing of the input data and output results for 
the NASA Task-1 transonic compressor with test conditions 
noted. TO compare the accuracy of the predicted flow with 
actual laboratory observations, computed axial velocit 
profiles at the rotor inlet, rotor outlet, stator inlet, and 
Stator outlet were compared with experimental results. in 
addition numerical results from Ref.7 were also compared. 

Pig 13-16 show the computer predictions plotted with the 
experimental values and the numerical solutions obtained by 
Hirsch and Warzee. The profiles shown were obtained after 
ten iterations and using a relaxation factor of 0.2. The 
figures show that the best overall agreenent with 
experimental data occurred in the stator inlet and outlet. 
In this region the worst error was 17% which occurred at the 
Stator tip inlet. The average error throughout the stator 


region with respect to experimental data was 6.6%. 


The rotor hub and tip outlet area exhibited 
instabilities in density convergence using equation 
PLease aie os Specifically, the density solution converged to 
within 8% at the rotor outlet tip and hub. It was found 
that by not allowing the nodal density at these nodes to go 
below a critical value of 0.06 lbm/cu ft, the solution for 
the stream function converged. By allowing the nodal 
densities at the rotor outlet tip and hub to go below this 
critical value, the computed velocities at these nodes 
became increasingly large and the arguement within the 
brackets of equation III.C.5.9 became less than one. This 
prevented continuation of the iterations for the stream 
function solution. In addition, the r©oter tip outlet 
exhibited more instability than the rotor hub outlet. The 
static density at the rotor hub outlet oscillated about a 
value of 0.062 lom/cu ft while the rotor tip outlet was 
constantly driven to the critical value of 0.06 lbm/cu ft. 
One method attempted to alleviate this problena was the 
following. Since a half-interval iteration routine was 
used, one trial run involved reversing the direction of 
consecutive guesses when the density iteration did not 
converge. It was found however, that after threa to four 
iterations of the systen of equations, the static densities 
at the rotor outlet tip and hub were again driven to smaller 
and smaller values which led to instability once more. The 


ne 


ct 


nodal densities converged at all interior points of 
rotor edge and mid-blade regions and also at all the rotor 
inlet nodes. By including all rotor nodes, the average 


error with respect to experimental data was 27.5%. 


Fig 17 shows a plot of convergence criteria,€, versus 
the number of iterations for a relaxation factor of 0.2. 
The stability of convergence is shown to initially decrease 
and then after the third iteration oscillates about an 
approximate value of 28%. It is important to nota that this 


curve represents the maximum value of € as shown in equation 
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(II.C.2.5). In addition, the curve in actuality represents 
the oscillation of nodal stream function values in the 
rotor/stator regions since in fact this is where the 


non-linearity is the greatest. 
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V. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STODY 


Agreement with both experimental data and numerical 
solution of Ref.7 was best in the stator region to within 
8%. Predicted axial velocity profiles in the rotor inlet 
area were within 26.2% of experimental results. The 
instabilities with respect to static density solutions are 
prevalent. One of the reasons for this numerical 
disagreement with Hirsch and wWarz2ee is the isentropic 
assumption imposed by the present program. Reconumendations 
for further study on the project include the addition of 
entropy variations in the rotor and stator blade regions. 
This would necessitate the use of blade correlation data 
{Ref.13] for loss predictions and involve additional input 


data plus program additions to Subroutine’s SLINE a'nd FCAL. 
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Figure 13 - AXIAL PROFILE AT ROTOR INLET 
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Figure 14 - AXIAL PROFILE AT ROTOR OUTLET 
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APPENDIX A 
COMPUTER PROGRAM 
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APPENDIX D 


CALCULATION OF ROTOR ELEMENT FLOW ANGLES 


The following is a brief synopsis of the procedure 
contained in Ret.13 for calculating the outlet relative 
flow angles in a rotor 2lement from the given inlet relative 
flow angle and blade solidity. The reader is referred to 
Ref.13, Chapter VI, for specific details of low speed 
correlation data. 


As stated in Section III.A, uniform flow conditions at 
the rotor blade edges were assumed. This assumption coupled 
with knowledge of the mass flow rate and rotational speed, 
enables one to calculate the inlet relative flow angl2.B,, 
as shown in Fig 18. 


Prom blade geometry information, the blade solidity,T, 
T= s (1) 


is obtained. At this point, ae and € are given and one may 
calculate Po , the rotor outlet relative flow angle from 
correlation curves depicted in Ref 13. The equation used to 
determine Pa is the following, 


Ba= Ke +8 (2) 


where K, is the angle between the tangent to the blade mean 


2 


117 


camber line and the axial direction (Fig 18). This is 
obtained from the blade geometry data. is the lew speed 
deviation angle which is obtained from the correlation 
curves in Ref. 13. The following equations show the 
relationship between S and the correlation data. 


f= §, +m (3) 


. > (ke). (KS), (a (4) 


The variables a, KS), aS), and ‘i. , are all values 
which are obtained from the correlation curves and are all 
functions of the given blad2 geometry. The quantity, > y is 
the blade camber angle and again is obtained from the blade 
geometry data. Once all the variables are obtained from the 
correlation data, equation (4) is solved for the deviation 
angle for an uncambered blade section, §., and then equation 
(3) is solved for the deviation angle, ¢ . One now 
calculates B, from equation (2) for the blade element. With 

Ba now a known quantity, one now calculates the absolute 
flow angle, d4,, from uniform flow assumptions. 


An example follows for node numbers 43 and 57 (Fig 46). 
From Ref. [12], Table II, the following quantities are 
obtained assuming the angle of incidence, i, (Fig 18) is 
zero and therefore the inlet relative flow angle, Bi , AS 
equal to k,. 


B = k= 61.88" 
T = 1.3062 
= 6.95° 


1138 


X, = 54.93° 
=) sax = 0.035 


tip radius = 18.25 in 


hub radius ot2S £0 
Assuming uniform flow at the rotor inlet and a 

rotational speed of 4359.5 RPM, the following quantities are 
determined from the rotor inlet velocity diagram (Fig 198). 

(107. & lbw /sec) (144 in*/FE) 

(0.08 \pm/et?) 1 (18.257 - 9.125") in? 

Vm = 246.802 Ft/sec 

where the area A, is determined from the hub and tip radii 


mM 
— 


Vo - ee, 


and the density is assumed to be 0.08 lbm/cu ft. 


Now one is ready to obtain the correlation data. From 
Ref.[13], Pig 162, with P, = 61.88° and T = 1.3062, 
oe 
§.) = 2650 


From Ref.(13], Pig 162, with & = 61.38° and © =1.3062, 


m= 0.235 
From Ref.{13], Fig 172, with t/c)max = 0.0350, 


K§)t = 0.29 
From Ref.{13], page 222, one uses the following value of 
(K§) sh for 65-series blades, 


K§)sh = 1.0 
At this point all the necessary data has been obtained for 
equations (3) and (4), 


a ° 
§. = (1.0) (0.29) (2.50) = 0.725 


°o 
From equation (3), 


§ = 0.725°+ 0.235(6.98) = 2.36° 
Finally, equation (2) gives the desired value of Bs + 


B, = 54.93 + 2.36 = 57.29° 
At this point the relative flow angle for node 57 has been 
obtained, = 57.29°. These two values of relative flow 
angles, f, = 61.88° for node 43 and 8, = 57.29" for node 57, 
are then read in the program as input data for numerical 


computation. 


This process is repeated at each required blade element 


section for the proper outlet relative flow angle. 
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Figure 18 - NOMENCLATURE FOR CASCADE BLADE 
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